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Abstract 

Concrete subjected to combined compressive stresses and temperature loading exhibits 
compressive strains, which are considerably greater than for concrete subjected to com- 
pressive stresses alone. This phenomenon is called transient thermal creep or load induced 
thermal strain and is usually modeled by macroscopic phenomenological constitutive laws 
which have only limited predictive capabilities. In the present study a meso-scale model- 
ing approach is proposed in which the macroscopically observed transient thermal creep 
results from the mismatch of thermal expansions of the meso-scale constituents. The 
meso-structure of concrete is idealized as a two-dimensional three phase material con- 
sisting of aggregates, matrix and interfacial transition zones (ITZ). The nonlinear ma- 
terial response of the phases are described by a damage-plasticity interface model. The 
meso-scale approach was applied to analyze compressed concrete specimens subjected to 
uniform temperature histories and the analysis results were compared to experimental 
results reported in the literature. 



Introduction 

The mechanical behaviour of concrete is significantly altered when exposed to high tem- 
peratures, whereby properties such as strength and stiffness are generally found to decrease 
with temperature and time dependent creep characteristics are found to be significantly 
enhanced. Experimental evidence also suggests that strength, stiffness and creep depend 
on the combined mechanical and thermal loading history. 
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When an unloaded concrete specimen is subjected to transient high temperatures it un- 
dergoes thermal expansion. When the same test is performed under sustained mechanical 
load, the total strains, minus the elastic strains, are significantly different from those ob- 
served in the unloaded case ( [Khoury, Grainger, and Sullivan 1985} IThelandersson 19871 
ISchneider IQSSP . This difference in the strains is normally referred to as transient thermal 
creep or load induced thermal strain. It is worth noting that the use of the word 'creep' 
is a subject for debate, since transient thermal creep is not time-dependent, as is the case 
for other forms of creep; for consistency with other work on this subject, the term 'creep' 
will continue to be used in this paper. 

Transient thermal creep is an irreversible thermo-mechanical strain, only observed on first 



heating and is considered a quasi-instantaneous response ( [Khoury, Grainger, and Sullivan 1985 
IThelandersson 1987^ . It is generally recognized that its presence provides an important 
relaxation of compressive stresses when concrete is heated beyond 100 °C. 

However, the irreversibility of transient thermal creep can lead to the build-up of poten- 



tially critical tensile strains in restrained concrete members during cooling (Nielsen, Pearce, and Bicanic 



Recent experimental studies by Cohna et al. (IColina and Sercombe 20041 IHassen and Colina 2006[ 



Sabeur, Meftah, Colina, and Platret 2008D performed at very slow heating rates up to 
temperatures of 400 °C have provided further confirmation of the existence of transient 
thermal creep as well as the principal characteristics outlined above. They also hypothe- 
sized that the phenomenon depends only on irrecoverable physical and chemical processes 
that take place due to increase of temperature. (IColina and Sercombe 2004jl also iden- 
tified that transient thermal creep is less in ordinary concrete (OC), which is character- 
ized by higher permeability, compared to that of high performance concrete (HFC) with 
lower permeability. (IHassen and Colina 2006t Sabeur, Meftah, Colina, and Platret 2008 ) 
concluded that dehydration of the cement paste is one of the critical processes in the 
development of transient thermal creep. 

There have been a number of efforts to phenomenologically model transient thermal creep 
strains during combined thermal and mechanical loading. Thelandersson (IThelandersson 1982p 
and Schneider (ISchneider 1988^ suggested uniaxial expressions for the transient thermal 
creep, with the former applying the theory of plasticity and extending the model to 
3D. Khennane and Baker (IKhennane and Baker 1992| IKhennane and Baker 19931) also 
adopted the theory of plasticity and based their numerical implementation of transient 
thermal creep on the work of de Borst and Peeters (Ide Borst and Peeters 19891) . The- 
landersson (IThelandersson 1987^ proposed a general 3D model that has been adopted by 
a number of authors (e.g. ( Nechnech, Meftah, and Reynouard 2002 )), but it is frequently 
unclear how the critical transient thermal creep parameter has been determined. Nielsen 
et al. ([Nielsen, Pearce, and Bicanic 2002|) derived the transient thermal creep parame- 



ter from experimental data reported by Schneider (ISchneider 1988^ and Khoury et al. 



(Khoury, Grainger, and Sullivan 1985). This model was subsequently implemented into a 



continuum Finite Element formulation by Pearce et al. (Pearce, Nielsen, and Bicanic 2002). 



Gawinetal. ( Gawin, Pesavento, and Schrefier 2004 Gawin, Pesavento, and Schrefier 2006 ) 
adopted a similar approach, although they related the transient thermal creep strains 



2 



to a thermo-chemical damage parameter rather than temperature. Bazant and Chern 
(IBazant and Chern 1987^ extended their original creep model to include the effect of tem- 



perature alterations, suggesting that both heating and cooling result in transient thermal 



creep. Thienel and Rostasy (Thienel and Rostasy 1996) proposed a general constitu- 



tive model for transient thermal strains based on their own experimental biaxial results. 
Youssef and Moftah (lYoussef and Moftah 20071) have undertaken a comparison of a num- 
ber of phenomenological models for transient thermal creep. 

In general, the mechanics of concrete subjected to temperature loading are not well un- 
derstood and the macroscopic phenomenological constitutive models which are based on 
existing experimental results cannot be relied upon to predict the structural response 
outside the range of these experiments. Furthermore, there is still dispute about the 
importance of moisture, drying shrinkage and differences in the thermal properties of 
the constituents of the meso-structure. The authors believe that a detailed study of the 
meso- and micro-structural behavior of concrete can help in resolving this matter by al- 
lowing the mechanisms which dominate the macroscopic response to be isolated. Recent 
promising research in this direction explored the macroscopic transient thermal creep of 
concrete via a two-phase meso-scale modeling approach which considered the thermal ex- 



pansion and the thermal damage of the elastic properties (Willam, Rhee, and Shing 2004 



Willam, Rhee, and Xi 2005 Grondin, Dumontet, Ben Hamida, Mounajed, and Boussa 2007) 



In the present work, a different meso-scale description is adopted, in which the meso- 
structure is enhanced by an additional phase, which represents the interfacial transition 
zone (ITZ) between the aggregates and the matrix. The macroscopic transient thermal 
creep is described by the nonlinear mechanical response, which is caused by the mismatch 
of thermal expansion of the different phases. The aim of the present study is to show that 
this mismatch of expansion has an important influence on the macroscopically observed 
transient thermal creep. 

The long term aim of the present research is to improve the understanding of the macro- 
scopic response of concrete subjected to temperature loading by investigating the influence 
of different processes on the meso and micro-scale. The present manuscript is focused on 
the nonlinear response due to the mismatch of thermal expansion on the meso-scale. 
Further work will be required to investigate other influences, e.g. moisture, creep and 
relaxation on micro and meso-scale. The present study is only a first step towards an 
understanding of the mechanisms that cause transient thermal creep. 

Meso-scale modeling of strongly heterogeneous materials like concrete is computationally 
very demanding. Continuum formulations for the description of fracture, especially at 
the interface between materials of significantly different stiffness, is difficult and can be 
accompanied by severe numerical problems. Therefore, a lattice approach is applied in the 
present study (IBolander and Saito 1998|IKawai 1977t Morikawa, Sawamota, and Kobayashi 1993 ) 



The heterogeneity of the material is described by spatially varying material proper- 
ties of the lattice elements with respect to their position within the meso-structure 
(Schlangen and van Mier 1992 ILilliu and van Mier 2003]) . which is composed of aggre- 
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gates, matrix and ITZs. It is recognised that each of the three phases comprise a com- 
plex microstructure. However, in the present study the response of these phases is de- 
scribed by phenomenological constitutive models. Alternatively, for at least some of 



these phases, micromechanics approaches could be used (Budiansky and O'Connell 1976 



Heukamp, Lemarchand, and Ulm 2005 Dormieux, Kondo, and Ulm 2006). 



Modeling approach 



The present approach to modeling the response of concrete subjected to combined ther- 
mal and mechanical loading is based on a meso-scale description. Aggregates, matrix 
and interfacial transition zones are modeled as separate phases with different material 
properties. A lattice approach is used to discretize the domain. For each lattice element, 
a nonlinear stress-strain relationship based on a combination of damage and plasticity is 
used. This damage-plasticity constitutive model was originally proposed for the modeling 



of concrete subjected to cyclic loading (Grassl and Rempling 2008) and is extended in the 



present study to describe the influence of thermal expansion. 



Lattice framework 



The domain is decomposed into polygons based on the Voronoi tessellation (lAurenhammer 199 ip 
The nodes for the Voronoi tessellation are placed sequentially in the domain, whereby the 
coordinates of each node are determined randomly. A minimum distance (im is enforced 
iteratively between the nodes. For each randomly placed node the distance to the exist- 
ing nodes is checked. Only if the smallest distance between two nodes is greater than the 
prescribed minimum distance, is the node accepted (IZubelewicz and Bazant 19870 . 

For this iterative process, the relationship between the number of nodes for a chosen 
domain and the minimum distance determines the distribution of the distances between 
nodes. The number of nodes for a minimum distance and domain size is expressed in 
the form of a density 

P=^ (1) 

If the density of the nodes is small, the variation of Voronoi cell sizes is large. On the 
other hand, if the domain is saturated with nodes, i.e. maximum number of nodes for 
a specific domain and minimum distance, the Voronoi cells are of similar size. Bolander 
and Saito (IBolander and Saito 1998^ determined numerically the density for a saturated 
node arrangement to be p ~ 0.68. The dual to the Voronoi tessellation is the Delaunay 
triangulation shown in Figure [1^. The edges of the Delaunay triangulation are the lattice 
elements, which describe the interaction between the nodes. 
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(a) 



(b) 



Figure 1: (a) Voronoi tessellation and dual Delaunay triangulation. Solid lines represent 
the cells defined by the Voronoi tessellation. Dashed lines represent the lattice elements 
describing the interaction between the nodes, (b) Meso-structure of concrete modeled by 
placing lattice elements perpendicular to the interfaces between aggregates and matrix. 



The use of an irregular lattice for the description of the fracture process of quasi-brittle ma- 
terials is important, since the fracture patterns are strongly influenced by the alignment of 
the lattice (IJirasek and Bazant 19951 Schlangen and Garboczi 1996). Related problems 



are known from continuum fracture approaches, where mesh alignment can have a strong 
influence on the fracture patterns obtained (IGrassl and Jirasek 2004t Grassl and Rempling 2007 
IJirasek and Grassl 2008]) . 

The meso-structure of concrete is discretized by placing lattice elements perpendicular to 
the interface between aggregates and the matrix. In Figure[T]3, a schematic meso-structure 
and the corresponding discretization is shown. 

The distribution of the aggregates is determined randomly using a cumulative distribu- 



tion function (Carpinteri, Cornetti, and Puzzi 2004). A standard pseudo-random number 



generator is used to generate probabilities from which the diameter d is determined. This 
procedure is repeated until the chosen volume fraction of aggregates is obtained. In 
the present 2D idealisation, aggregates are idealized as cylindrical inclusions of diameter 
d. The aggregates are placed sequentially by means of randomly generated coordinates 
within the area of the specimen. For each set of generated coordinates it is checked 
that no overlap with existing aggregates occurs. However, overlaps with boundaries are 
permitted. 

The lattice elements describe the elastic and inelastic response of concrete, i.e. both 
the elastic and inelastic displacements are described by the displacement jump between 
rigid bodies. Each rigid body possesses two translations and one rotation. In the 
local coordinate system (n and s), shown in Figure O the degrees of freedom Ue = 
{ui,vi, (pi,U2,V2, 4>2}'^ of two rigid bodies sharing an interface describe the displacement 
discontinuities Uc = {uc, Vc}^ at the mid point C of the interface, where the springs are 
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Figure 2: a) Lattice element representing the connection between two nodes, b) Geometry 
in the local co-ordinate system. 

element axis and h is the length of the element (Figure [2]). If C is on the left hand side 
of the element, e is positive. Otherwise, e is negative. The displacements Uc at the point 
C are transformed into strains £ = {e^i 6^}'^ = Uc/h. The strain £ is related to the stress 
(T = {(Tn cTg}^ by the constitutive model of the springs, which is described in the following 
section. 

The stiffness matrix of a lattice element in the local coordinate system results in 

Kt = ^B^DtB (4) 

n 

where A is the cross sectional area of the common interface of the two rigid bodies. 



Constitutive model 



In the present study the response of concrete is analyzed with an interface model, which 
was developed in ( Grassl and Rempling 2008 ) for cyclic loading and extended here to 
include thermal loading. In the following, the modeling approach is briefly reviewed to 
introduce the parameters, which are investigated in the present study. The interface model 
relies on stress based elasto-plasticity. The strains are related to the stress cr = (an, cts)^ 
as 



(5) 
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where De is the elastic stiffness, Sp = (£:pn,£^ps)^ is the plastic strain. The mechanical 
strain is defined as 

£m = £ - £t (6) 

where the thermal strains are 

St = aT (T - (^^ (7) 

The parameter T is the temperature within the specimen and Ta is the reference temper- 
ature. The thermal expansion coefficient ot is chosen as a function of the temperature: 

ot (T) = Oti + «T2(r - Ta) + OTsiT - T^f (8) 

where oti, aT2 and axs are model parameters. 

The elastic stiffness is 

D= - { 4} («' 

where E and 7 are model parameters controlling both the Young's modulus and Poisson's 
ratio of the material (IGriffiths and Mustoe 200TI) . 

The small strain plasticity model consists of the yield surface, flow rule, evolution law for 
the hardening parameter and loading and unloading conditions. A detailed description 



of the components of the plasticity model is presented in ( [Grassl and Rempling 20()8| ). 
The initial yield surface is determined by the tensile strength /t, by the ratio s of the 
shear and tensile strength, and the ratio c of the compressive and tensile strength. The 
evolution of the yield surface during softening is controlled so that linear stress inelastic 
displacement laws for pure tension and compression are obtained, which are characterized 
by the fracture energies Gft and The eight model parameters 7, /t, s, c, Gft and Gfc 
can be determined from a tensile, shear and compressive test of the material phase. The 
constitutive response of the interface model is demonstrated by the stress-strain response 
for fluctuating normal strains in Figure [31 The normal strain is increased to point A. 
Then the strain is reduced to point B and again increased to point C. 



Structural analysis 

The lattice modeling approach is applied to the meso-scale analysis of the macroscopic 
phenomenon of transient thermal creep. Experimental results performed by Thelander- 
sson and reported in (IThelandersson 1987^ show that there is no significant influence of 
the rate of temperature increase for rates of l°C/min and 5°C/min. This indicates that 
the temperature loading in these experiments is slow enough to neglect the influence of 
the thermal capacity of the materials. Therefore, the thermal analysis was simplified 
by assuming a uniform temperature increase throughout the specimen, disregarding the 
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normalised strain (e^ ft)/(2 G(t) 



Figure 3: Stress-strain response for fluctuating normal strains for /i = 1 (solid line) and 
/i = (dashed line). 

thermal conductivity and capacity, which may differ significantly for the three phases and 
may play an important role for increased heating rates. 

A concrete prism was idealized by a two-dimensional discretization under the assumption 
of plane stress. The height and width of the specimen was chosen as L = 300 mm and 
D = 150 mm, respectively. The mesh was generated with a vertex density of p = 1.5 
and a minimum distance of = 3 mm. The distribution of aggregates distribution was 
obtained with c/max = 32 mm, (imin = 10 mm and an aggregate volume fraction of pa = 0.3. 
Aggregates of a diameter less than 10 mm are not modelled discretely. The geometry, 
loading setup and mesh is shown in Figure HI 

Three types of loading situations were investigated. Firstly, the concrete specimen was 
subjected to monotonic uniaxial compression at ambient temperature. Secondly, the 
specimen was heated while subject to various levels of constant uniaxial compressive 
stress. Thirdly, the specimen was subjected to uniform heating under uniaxial restraint. 
The results of the latter two loading scenarios were compared to experimental results 
obtained by Thelandersson and reported in (IThelandersson 1987p . 

Uniaxial compression 

The uniaxial compression analysis was controlled by the top displacement Am. The aver- 
age strain is determined as e = Au/L. The average stress a is determined as P/D, where 
P is the reaction force caused by the top displacement Au. 
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Figure 4: Model of a concrete prism: a) Geometry and loading setup, b) Delaunay tri- 
angulation c) Voronoi tessellation d) Interfaces representing the ITZ between matrix and 
aggregates. 



Table 1: Model parameters for the structural analysis. 



E [GPa] 7 ft [MPa] q c Gft [J/m^] Gfc [J/m^] /i 

Matrix 24 0.33 6 2 20 600 600000 1 

ITZ 34.3 0.2 1.5 2 20 150 150000 1 

Aggregates 60 0.08 _ _ _ . . _ 



The model parameters for the mechanical part of the model are presented in Table [H 
The parameters were chosen to obtain a macroscopic elastic stiffness of = 21.5 GPa, 
a compressive strength of fc = 35 MPa and an inelastic strain at peak of ^i^max = 0.0015, 
which were reported in (IThelandersson 1982^ . The stress-strain response is shown in 
Figure [5l The fracture patterns obtained at four stages of the analysis marked in Figure [5] 
are shown in Figure El 

At an early stage within the pre-peak regime, the cracks are mainly restricted to the 
aggregate-matrix interface and are perpendicular to the loading direction (Figure [6^). 
Around the peak load (Figures and c), the interfaces between aggregates and matrix 
become strongly damaged and cracks occur in the matrix phase, which are orientated 
parallel to the direction of the compressive stress. A similar direction of crack-propagation 



was observed in (Pichler, Hellmich, and Dormieux 2007) for the case of penny-shaped 



cracks embedded in an elastic matrix. At a later stage the cracks localize into a shear 
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Figure 5: Stress-strain response in uniaxial compression. 




Figure 6: The fracture patterns at four stages of analyses marked in Figure O The black 
lines represent interfaces for which the damage variable increases. 



The plasticity constitutive model limits the strength in tension, shear and compression 
(Figure [3]). In the uniaxial compressive test, failure in the matrix is entirely caused by 
tensile/low shear stress states. However, the compressive strength in the ITZs is reached 
in the pre-peak regime of the macroscopic load-displacement curve. 
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Table 2: Thermal expansion coefficients for the structural analysis. 

axi [°C-^] aT2 [°C-^] aT3 ["C"^] 
Matrix 1.5 x 10"^ -6 x 10"^^ 

ITZ 1.5 X 10"^ -6 X 10-11 

Aggregates 1.5 x 10~^ 7 x 10"* 



Heating under various compressive loads 

The second group of tests consists of two steps. First a uniaxial compression stress is 
applied. Then, this stress is kept constant and the temperature is increased from 20 °C 
to 1000 °C, while the strain in the axial direction is recorded. These tests are carried out 
for four stress levels of a/ fc = 0, 0.225, 0.45 and 0.675. 

The thermal expansion coefficients for the matrix, interface and aggregates are presented 
in Table [2J The Matrix and ITZ phase initially expands and then contracts for higher 



temperatures (Khoury, Grainger, and Sullivan 1985), corresponding to dehydration of the 



phases. All other parameters are the same as in the previous section. The thermal 
expansion coefficients of the three phases were determined from a ffi of the model response 
to the experimental results for free expansion (cr//c = 0). All other analyses are carried 
out with these parameters and no additional fits were performed. Thus, the analysis 
for cr//c = 0.225, 0.45 and 0.675 are model predictions. The strain-temperature curves 
obtained with the meso-scale approach are compared to the experimental curves reported 
in (IThelandersson 1987^ in Figure [3 The fracture patterns for the four different loading 



cases at a late stage of analysis marked in Figure [7] are shown in Figure El 

For the case of free expansion (o"//c = 0) cracks in the matrix propagate radially outwards 
from the aggregate inclusions due to the thermal incompatibility of the different phases. 
For the other analyses, the crack orientations are influenced by the applied compressive 
stresses. For a stress level of cr//c = 0.675 the majority of the damaged interfaces in the 
matrix phase are parallel to the loading direction. 



Heating under uniaxial restraint 

The specimen is restrained in the axial direction and the temperature is increased from the 
ambient temperature 20 °C to 1000 °C. The same material properties as in the previous 
example were used. The normalized average stress obtained from the meso-scale analysis 
is plotted against temperature in Figure [9] and compared with the experimental results 
of Thelandersson reported in (IThelandersson 1987^ . Initially, the compressive stresses 
increase with increasing temperature until approximately 200 °C and then reduce with 
further temperature increase. It should be noted that the model parameters adopted 
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Figure 7: Strain versus temperature for constant compressive stresses (a) a/ fc = 0, (b) 
a//c = 0.225, (c) a//c = 0.45 and (d) a/f^ = 0.675. 




Figure 8: The fracture patterns for (a) a/ fc = 0, (b) 0.225, (c) 0.45 and (d) 0.675 at late 
stages of analyses marked in Figure [3 The black lines represent interfaces for which the 
damage variable increases. 

are not fitted to match the experimental results. Nevertheless, the meso-scale modeling 
approach in this work qualitatively captures the response observed in the experiments. 
The fracture patterns at two stages of the analysis marked in Figure [9] are shown in 
Figure [TDl 
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Figure 9: Stress versus temperature for heating under uniaxial restraint. 




Figure 10: The fracture patterns at two stages of analysis for heating under uniaxial 
restraint marked in Figure M The black lines represent interfaces for which the damage 
variable increases. 

Macroscopic thermal damage 

Finally, the capabilities of the meso-scale approach are further illustrated by evaluating 
the stiffness at different stages of heating without any restraints. The temperature is 
increased and a small uniaxial stress increment is applied at different temperature levels 
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to determine the stiffness of the material in the axial direction. The macroscopic stiffness 
-Ejn is computed as 

Aa , , 

= ^ (10) 

The stiffness obtained at ambient temperature is used as reference stiffness i?o to compute 
the macroscopic thermal damage during heating as 

Eq — Ej^ 

Ut = (11) 

The evolution of the macroscopic thermal damage during heating is presented in Figure fTTl 




200 400 600 800 1000 

temperature T [°C] 



Figure 11: Thermal damage versus temperature. 



Conclusions 



This paper presents a three-phase meso-scale approach to modelling concrete subjected 
to thermo-mechanical loading. An elasto-plasticity model for the mortar and interfacial 
transition zones is adopted that is able to describe the fracture process of concrete in com- 
pression. The model is also capable of qualitatively describing the main characteristics of 
concrete subjected to combined thermo-mechanical loading and the results indicate that 
the mismatch of thermal expansion has a strong influence on the thermal transient creep. 
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Moreover, the path-dependence of the thermo-mechanical loading is described quahta- 
tively correctly by the model. It is worth reiterating that, unlike other approaches, no 
phenomenological thermal damage mechanism is introduced to model the experimentally 
observed results, although phenomenological models for mechanical damage and thermal 
expansion of the individual phases are required. The present study is only a first step in 
a longer research project and further studies are required to investigate the influences of 
moisture and basic creep on thermal transient creep and the results will be presented in 
future publications. 
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